/*    Copyright 2010 Tobias Marschall
 *
 *    This file is part of MoSDi.
 *
 *    MoSDi is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 3 of the License, or
 *    (at your option) any later version.
 *
 *    MoSDi is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with MoSDi.  If not, see <http://www.gnu.org/licenses/>.
 */

package mosdi.tests;

import junit.framework.TestCase;

import mosdi.distributions.PoissonDistribution;

import org.junit.Test;

public class PoissonDistributionTest extends TestCase {

	@Test
	public void testGetProb() {
		PoissonDistribution pd = new PoissonDistribution(4.5);
		assertEquals(1.0, pd.getProb(0)/1.1108996538242306e-02d, 1e-10);
		assertEquals(1.0, pd.getProb(17)/3.9739193639645268e-06d, 1e-10); 
		assertEquals(1.0, pd.getProb(100)/2.4941499169512807e-95d, 1e-10); 
	}

	@Test
	public void testGetQuantileByPValue() {
		PoissonDistribution pd = new PoissonDistribution(3.35);
		double threshold = 1e-85d;
		int k = pd.getQuantileByPValue(threshold);
		double d = 0.0;
		for (int i=k; pd.getProb(i)>0.0; ++i) d+=pd.getProb(i);
		assertTrue(d>threshold);
		d = 0.0;
		for (int i=k+1; pd.getProb(i)>0.0; ++i) d+=pd.getProb(i);
		assertFalse(d>threshold);
	}

	@Test
	public void testCCDF() {
		PoissonDistribution pd = new PoissonDistribution(3.5);
		double[] expected = {1.0, 9.698026e-01,8.641118e-01,6.791528e-01,4.633673e-01,2.745550e-01,1.423864e-01,6.528810e-02,2.673892e-02,9.873658e-03,3.314944e-03,1.019394e-03,2.889922e-04,7.595824e-05,1.860294e-05,4.264112e-06,9.183859e-07,1.865083e-07,3.582763e-08,6.528614e-09,1.131426e-09,1.869184e-10,2.950038e-11,4.456613e-12,6.456050e-13,8.983291e-14,1.202482e-14,1.550659e-15,1.928963e-16,2.317603e-17,2.692547e-18,3.028075e-19,3.299816e-20,3.487763e-21,3.578728e-22,3.567819e-23,3.458735e-24,3.262881e-25,2.997549e-26,2.683543e-27,2.342660e-28,1.995428e-29,1.659371e-30,1.347953e-31,1.070192e-32,8.308567e-34,6.310737e-35,4.691661e-36,3.415547e-37,2.435946e-38,1.702661e-39};
		double[] actual = pd.getCCDF(51);
		for (int i=0; i<=50; ++i) {
			assertEquals(1.0, expected[i]/actual[i],1e-6);
		}
	}

//	@Test
//	public void testGetQuantileByPValue2() {
//		Log.getInstance().setTimingActive(true);
//		Log.getInstance().startTimer();
//		FactorialCalculator fc = new FactorialCalculator(1000);
//		for (int i=0; i<10000; ++i) {
//			double lambda = 0.005*i;
//			PoissonDistribution pd = new PoissonDistribution(lambda, fc);
//			double threshold = 1e-50d;
//			int k = pd.getQuantileByPValue(threshold);
////			System.out.println(String.format("lambda=%f ==> k=%d", lambda, k));
//		}
//		Log.getInstance().stopTimer("zzz");
//	}

	@Test
	public void testCompoundDistribution() {
		PoissonDistribution pd = new PoissonDistribution(5.0);
		double[] clumpSizeDistribution = {0.0, 0.0, 0.9, 0.0, 0.1};
		double[] actual = pd.compoundPoissonDistribution(clumpSizeDistribution, 99);
		assertEquals(100, actual.length);
		double[] expected = {0.006737946999085467, 0.0, 0.030320761495884602, 0.0, 0.071590686865283082, 0.0, 0.11749295079655281, 0.0, 0.15007724136244271, 0.0, 0.15856810738550903, 0.0, 0.14393895409953877, 0.0, 0.11518477154763336, 0.0, 0.082783803257986122, 0.0, 0.054190209578730067, 0.0, 0.032663974636227118, 0.0, 0.018288917767431995, 0.0, 0.009580342049139259, 0.0, 0.0047231120760429721, 0.0, 0.0022024533136666157, 0.0, 0.00097561013250284953, 0.0, 0.0004120436818705899, 0.0, 0.00016645921770120617, 0.0, 6.4506120084778745e-05, 0.0, 2.4038776741195285e-05, 0.0, 8.6340307710078772e-06, 0.0, 2.9948531052728932e-06, 0.0, 1.0050395338516324e-06, 0.0, 3.2684917424370616e-07, 0.0, 1.0316086741451294e-07, 0.0, 3.1642923104360566e-08, 0.0, 9.4443854378513568e-09, 0.0, 2.7460243546182074e-09, 0.0, 7.7862482262975959e-10, 0.0, 2.1551158815352131e-10, 0.0, 5.8280898977353449e-11, 0.0, 1.54121172113423e-11, 0.0, 3.9886070758873009e-12, 0.0, 1.0109348197828828e-12, 0.0, 2.5111216955618436e-13, 0.0, 6.1169702365306058e-14, 0.0, 1.4621550838890603e-14, 0.0, 3.4315319227111824e-15, 0.0, 7.9114327608133964e-16, 0.0, 1.7927376064300527e-16, 0.0, 3.994687997437155e-17, 0.0, 8.7569444031140687e-18, 0.0, 1.8893602330567803e-18, 0.0, 4.0137361515975717e-19, 0.0, 8.3989579574447345e-20, 0.0, 1.7318371627661534e-20, 0.0, 3.520048954324436e-21, 0.0, 7.055019557898183e-22, 0.0, 1.3947516157038781e-22, 0.0, 2.7206942507276801e-23,6.46308397646e-24};
		for (int i=0; i<100; ++i) {
			if (expected[i]==0.0) assertTrue(expected[i]==0.0);
			else assertEquals(1.0, expected[i]/actual[i],1e-10);
		}
	}

	@Test
	public void testCompoundLogDistribution() {
		PoissonDistribution pd = new PoissonDistribution(5.0);
		double[] clumpSizeDistribution = {0.0, 0.0, 0.9, 0.0, 0.1};
		double[] actual = pd.logCompoundPoissonDistribution(clumpSizeDistribution, 99);
		assertEquals(100, actual.length);
		double[] expected = {0.006737946999085467, 0.0, 0.030320761495884602, 0.0, 0.071590686865283082, 0.0, 0.11749295079655281, 0.0, 0.15007724136244271, 0.0, 0.15856810738550903, 0.0, 0.14393895409953877, 0.0, 0.11518477154763336, 0.0, 0.082783803257986122, 0.0, 0.054190209578730067, 0.0, 0.032663974636227118, 0.0, 0.018288917767431995, 0.0, 0.009580342049139259, 0.0, 0.0047231120760429721, 0.0, 0.0022024533136666157, 0.0, 0.00097561013250284953, 0.0, 0.0004120436818705899, 0.0, 0.00016645921770120617, 0.0, 6.4506120084778745e-05, 0.0, 2.4038776741195285e-05, 0.0, 8.6340307710078772e-06, 0.0, 2.9948531052728932e-06, 0.0, 1.0050395338516324e-06, 0.0, 3.2684917424370616e-07, 0.0, 1.0316086741451294e-07, 0.0, 3.1642923104360566e-08, 0.0, 9.4443854378513568e-09, 0.0, 2.7460243546182074e-09, 0.0, 7.7862482262975959e-10, 0.0, 2.1551158815352131e-10, 0.0, 5.8280898977353449e-11, 0.0, 1.54121172113423e-11, 0.0, 3.9886070758873009e-12, 0.0, 1.0109348197828828e-12, 0.0, 2.5111216955618436e-13, 0.0, 6.1169702365306058e-14, 0.0, 1.4621550838890603e-14, 0.0, 3.4315319227111824e-15, 0.0, 7.9114327608133964e-16, 0.0, 1.7927376064300527e-16, 0.0, 3.994687997437155e-17, 0.0, 8.7569444031140687e-18, 0.0, 1.8893602330567803e-18, 0.0, 4.0137361515975717e-19, 0.0, 8.3989579574447345e-20, 0.0, 1.7318371627661534e-20, 0.0, 3.520048954324436e-21, 0.0, 7.055019557898183e-22, 0.0, 1.3947516157038781e-22, 0.0, 2.7206942507276801e-23,6.46308397646e-24};
		for (int i=0; i<100; ++i) {
			if (expected[i]==0.0) assertTrue(Double.isInfinite(actual[i]));
			else assertEquals(Math.log(expected[i]),actual[i],1e-10);
		}
	}

	@Test
	public void testCompoundPValue() {
		PoissonDistribution pd = new PoissonDistribution(10.3);
		double[] clumpSizeDistribution2 = {0.0, 0.0, 0.9, 0.0, 0.1};
		assertEquals(1.0, 0.567046831342/pd.compoundPoissonPValue(clumpSizeDistribution2, 21),1e-10);
		assertEquals(1.0, 6.5874302605e-174/pd.compoundPoissonPValue(clumpSizeDistribution2, 501),1e-10);
		pd = new PoissonDistribution(5.8);
		double[] clumpSizeDistribution = {0.0, 0.8, 0.1, 0.05, 0.0, 0.05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
		assertEquals(1.0, 8.47537652912e-64/pd.compoundPoissonPValue(clumpSizeDistribution, 200),1e-10);
	}

	@Test
	public void testCompoundLogPValue() {
		PoissonDistribution pd = new PoissonDistribution(10.3);
		double[] clumpSizeDistribution2 = {0.0, 0.0, 0.9, 0.0, 0.1};
		assertEquals(Math.log(0.567046831342), pd.logCompoundPoissonPValue(clumpSizeDistribution2, 21),1e-10);
		assertEquals(Math.log(6.5874302605e-174), pd.logCompoundPoissonPValue(clumpSizeDistribution2, 501),1e-10);
		pd = new PoissonDistribution(5.8);
		double[] clumpSizeDistribution = {0.0, 0.8, 0.1, 0.05, 0.0, 0.05, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
		assertEquals(Math.log(8.47537652912e-64),pd.logCompoundPoissonPValue(clumpSizeDistribution, 200),1e-10);
	}

	private static void assertPrecisionLimit(double expectation) {
		PoissonDistribution pd = new PoissonDistribution(expectation);
		int n = pd.pvaluePrecisionLimit();
		double[] clumpSizeDist = {0.0, 1.0};
		assertTrue(pd.compoundPoissonPValue(clumpSizeDist, n)<=Double.MIN_VALUE);
		assertTrue(pd.compoundPoissonPValue(clumpSizeDist, n-1)>=Double.MIN_VALUE);
	}
	
	public void testPrecisionLimit() {
		assertPrecisionLimit(20.0);
		assertPrecisionLimit(0.5);
		assertPrecisionLimit(17.5);
		assertPrecisionLimit(1024);
		assertPrecisionLimit(523.123);
		assertPrecisionLimit(777.777);
	}
	
}
